Deformation of a self-propelled domain 
in an excitable reaction-diffusion system 
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I. INTRODUCTION 



Takao Ohta,^ Takahiro Ohkuma,^ and Kyohci Shitara^ 
^ Department of Physics, Kyoto University, Kyoto 606-8502, Japan 

We formulate the theory for a self-propelled domain in an excitable reaction-diffusion system in 
two dimensions where the domain deforms from a circular shape when the propagation velocity is 
increased. In the singular limit where the width of the domain boundary is infinitesimally thin, 
' we derive a set of equations of motion for the center of gravity and two fundamental deformation 

' modes. The deformed shapes of a steadily propagating domain are obtained. The set of time- 

, evolution equations exhibits a bifurcation from a straight motion to a circular motion by changing 

the system parameters. 

< 

m 

@Self-organizcd dynamics of domains in reaction-diffusion media have attracted much attention recently [1-3]. 
' (-H ' Formation and propagation of domains in excitable reaction-diffusion systems have been studied both theoretically 
and experimentally. Not only extended domain patterns such as spiral wave but also interacting disconnected domains 
1^ , have been found in experiments of chemical reactions where fascinating dynamics due to propagation, collision and 
H ' self-replication of domains have been observed [4, 5]. Throughout the present paper, we shall use the word "domain" 
^ I , ' for a localized isolated object in two dimensions rather than the word "pulse" which is the concept in one dimension. 

Computer simulations of reaction-diffusion equations have revealed various interesting dynamics of domains [6-10]. 
^ Krisher and Mikhailov have investigated numerically collisions of a pair of domains in two dimensions in an excitable 
reaction-diffusion system with a global coupling [10]. They have also shown that an isolated domain is deformed 
' substantially from a circular shape when the propagating velocity is increased. It has been shown that Sierpinski 
gasket patterns appear in space-time plot of the pulse trajectory in one dimension in several reaction diffusion systems 
due to a delicate interplay between pair-annihilation and self-replication of pulses upon collision [11-13]. 

It has been known in excitable reaction-diffusion equations that a motionless localized domain loses its stability 
and begin to propagate when a system parameter exceeds the critical threshold. Reduction theories in the form of 
^ I ordinary differential equations in terms of a few degrees of freedom have been developed near this bifurcation called 
a drift bifurcation not only in one dimension [14, 15] but also in two and three dimensions [16-20]. For example, the 
motion of the center of gravity of an isolated spherical domain obeys in the vicinity of a supercritical drift bifurcation 

I> , — =71^- ji^l (1) 

0\ 

where v is the velocity. The coefficient 7 is negative below the threshold whereas it is positive above the threshold. 
The interaction between two domains can be incorporated in Eq. (1) by assuming that the characteristic length of 
the domain is much smaller than the distance between the domains. Here it is remarked that Eq. (1) for 7 > is a 
model equation for the active Brownian motion if a noise term is added [21-23]. 

As mentioned above, a domain is generally deformed when the propagating velocity is increased. However, to 
our knowledge, no analytical theory has been available for dynamics of deformable self-propelled domains. Only a 
stationary deformed shape of a steadily propagating domain at a constant speed has been studied by means of the 
interface dynamics in reaction diffusion systems [20, 24]. (In Ref. [20], the author claims that a stable propagating 
^ ' domain exists in a set of two-component reaction diffusion equations whereas the author in Ref. [24] studies a 
\ traveling domain in the slightly different reaction diffusion equations introducing a feedback control term to stabilize 
the traveling domain.) 

In the present paper, we shall derive the time-evolution equations for the motion of a deformed domain starting with 
the excitable reaction-diffusion equations studied numerically by Krisher and Mikhailov [10]. In a previous paper [25], 
we introduced phenomenologically a set of evolution equations for v and S, the latter of which is the so-called nematic 
order parameter tensor [26] to represent an elliptical deformation of domain. A remarkable property of this system 
is that there is a bifurcation such that a straight motion of a domain becomes unstable for sufhciently small values 
of the relaxation rate of S and the domain undergoes a circular motion. Therefore a deformable domain exhibits two 
bifurcations. One is the drift bifurcation and the other is the straight-circular-motion bifurcation. Hereafter, we shall 
call this bifurcation a rotation bifurcation. 

We shall employ the following assumptions in the derivation of the time-evolution equations for a deformable self 
propelled particle. First of all, we represent the domain dynamics by means of the interface equation of motion for 
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the domain boundary assuming that the interface width is much smaller than any other lengths of the problem [3]. 
The second assumption is slowness of the propagating velocity and hence weakness of deformation provided that the 
system is near the drift bifurcation point. 

Although the present formulation is specific to the excitable reaction-diffusion equations, the set of the time- 
evolution equations takes a general form in terms of the vector w, the second rank tensor S and the third rank tensor 
U . The tensor U was not considered in ref. [25] but is necessary to represent the head-tail asymmetry of a propagating 
domain. It is emphasized that the new tensor variable U does not enter up to the leading order in the equations for 
V and S and hence the previous analysis in ref. [25] is unaltered. 

We summarize some of the related studies. First of all, it should be remarked here that there are several experiments 
of spontaneously deformed propagating domains. For example, an oil droplet in surfactant solutions undergoes a self- 
propelled motion and changes its shape depending on the propagating velocity [27, 28]. The present theory would 
provide some insight into these experiments. Next, Shimoyama et al [29] have introduced a model of interacting 
elliptical particles, where the direction of the long axis is chosen to be an independent dynamical variable coupled 
to the velocity of the center of gravity. The reaction driven propulsion has been studied both experimentally and 
theoretically [27, 28, 30-34]. The active Brownian particles with a harmonic interaction between a pair of particles 
show a transition from a translational motion to a rotational motion by increasing the noise strength [35]. The 
Langevin equations for the center-of-mass position of rods and the rod orientation have been introduced to study 
the Brownian circular swimmer in a confined geometry [36] . Collective dynamics of self-propelled particles with an 
interaction between chemotactic materials and an internal degree of freedom has been investigated [37]. However, all 
of these theoretical studies deal with undeformable objects. 

The present paper is organized as follows. In the next section (section 2), we start with the reaction-diffusion 
equations for an activator and an inhibitor with a global coupling. In section 3, we briefly review the interface 
dynamics. The representation of a deformed domain is given in section 4. Apart from the center of gravity, we 
consider two tensor variables S and U corresponding, respectively, to the n — ±2, and ±3 modes in the Fourier series 
expansion of the deformation 5R{4>,t) = '-"(0^"'^" with (j) the angle to specify the location on the interface. In 
section 5, the time-evolution equation for the center of gravity is derived whereas, in section 6, the equations for 
n = ±2, and ±3 modes are obtained to complete the set of the equations for a deformable self-propelled domain. 
The approximations used in the derivation of the equations and the stationary shape of a propagating domain are 
described in section 7. Summary and discussion are given in section 8. 

II. EXCITABLE REACTION DIFFUSION SYSTEM 
We start with a coupled set of reaction-diffusion equations for an activator u and an inhibitor v. 

Te^ ^e^S/^u + f{u,v}~v , (2) 



^ = DV'^v + u--fv , (3) 

where 

f{u,v} = -u + e{u-p'{u,v}) , (4) 
with 9{x) = 1 for x > and 9{x) = for a; < 0. The functional p'{u, v} contains a global coupling as 

p' =p + (j[j(u + v)dr - W] , (5) 

where a and W are positive constants, < p < 1/2 and the integral runs over the whole space. The constants r and 
7 are positive and chosen such that the system is excitable and that a localized stable pulse (domain) solution exists. 
Inside the domain, the variable u is positive surrounded by the rest state where u and v vanishes asymptotically away 
from the domain. The parameter e is a measure of the width of the domain boundary. Hereafter we consider the limit 
0. 

The set of Eqs. (2) and (3) with (4) and (5) was introduced by Krisher and Mikhailov [10]. They investigated the 
collision of domains in two dimensions by computer simulations and found a reflection of a pair of colliding domains 

when the propagating velocity is sufficiently small. The reason why the global coupling is necessary in the system (2) 
and (3) is as follows. It has been known that a motionless domain is stable when r is sufficiently large. There is a 
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bifurcation such that below a threshold value a motionless domain becomes unstable and. when the global coupling 
is absent, a breathing motion occurs such that the radius of domain undergoes a periodic oscillation [38 40]. On the 
other hand, it is also known that the drift bifurcation exists for smaller values of r. Therefore, when the value of t 
is decreased, one generally encounters the situation that the breathing bifurcation occurs before the drift bifurcation. 
In order to study the intrinsic property of a propagating domain, one must exclude the breathing instability. This is 
achieved if one chooses a sufficiently large value of a in the global coupling (5) so that it becomes p' = p with 
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dr{u + v) = W , (6) 



and f{u,v} is no more a functional but is given by 

f{u) = -u + e{u-p) . (7) 
In the limit e and a — > oo in (5), Eq. (2) becomes 

-u + e{u-p) -V = . (8) 

The location of the domain boimdary (interface) is defined such that u = p and u > {<)p inside (outside) the domain. 
It is shown from Eq. (8) that u + 7; = 1 inside the domain whereas u + v — outside the domain. Therefore the 
constraint (6) means that the volume of an excited domain is independent of time and hence the breathing motion 
is prohibited. It is also noted that the bifurcation from a motionless domain to a propagating domain becomes 
supercritical in the limit a — > 00 [10]. 
Substituting Eq. (8) into Eq. (3) yields 

— = DV^v + 9{u-p)- f3v , (9) 

where /? = 1 + 7. The equilibrium solution of a motionless circular domain with radius Rq in two dimensions has been 
obtained in ref. [39]. Noting the relation that 9{u — p) = 9{Ro — r), the equilibrium solution w = w is given from Eq. 
(9) for < r < i?o by 

u = i[l - RoKKi{RoK)Io{rK)] , (10) 

and for r > i?o by 

V = ^hiRoK)Ko{rK) , (11) 

where 

^ (12) 

and /„ and are the modified Bessel functions. The equilibrium profile of u is given by the relation u = 0{Ro — r) — v. 

From Eqs. (10) and (11), it is clear that the inhibitor v changes gradually in space with the characteristic length 
K~^. On the other hand, the activator u is discontinuous at r = i?o in the limit e — > and hence there is a sharp 
interface. In this situation, one can derive the interface equation of motion as shown in the next section. 




III. INTERFACE EQUATION OF MOTION 

In order to see the spatial variation of u near r = Rq, one must rescale the space coordinate as r' = r/e ^ 0(1). In 

this length scale, the spatial variation of v is negligible and the value of v in (2) can be replaced by the value at the 
interface v{r,t) = w. As mentioned above, the location of the interface is defined through the condition u{r,t) = p. 
For a given value of w, one can readily obtain, from Eq. (2), the equation of motion for an arbitrary deformed interface 
[39] 



tV = eK + Tc{w) + L , 



(13) 
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where V is the normal component of the velocity directed from the inside to the outside of a domain and K is the 
mean curvature defined such that it is positive when the center of the curvature is outside the excited domain. The 
second term in (13) is the velocity for a flat interface and is related with w as 

CT 

l-2p-2w. (14) 



V(ct)2+4 



The unknown constant w is determined by solving Eq. (3) or Eq. (9) for a given interface configulation. The last 
term L is a Lagrange multiplier for the constraint of the domain area conservation (6). That is, the constant L is 
determined by the condition 
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dLoV = , (15) 



where doj is the infinitesimal length (area in three dimensions) of the interface and the integral runs over the interface. 

When the motion of the interface is slow compared with the relaxation rate of the inhibitor, one may deal with the 
term dv/dt in (9) as a perturbation so that the asymptotic solution of (9) can be written as 

where G is defined through the relation 

{-DV^ + l3)G{r - r') = S{r - r') , (17) 

and the abbreviation such that OA = J dr'G{r — r')A{r') has been used. The inhibitor v given by Eq. (16) at the 
interface position determines the value of w. Substituting it into Eq. (13) with (14) completes the closed form of the 
interface equation of motion. As discussed in ref. [18], the motion of interface is arbitrarily slow in the vicinity of the 
supercritical drift bifurcation where the expansion (16) is justified. 

IV. DEFORMED DOMAIN 

We consider a deformed domain with the center of gravity p{t). Its time-derivative p is given by 

p=v = ^j du)V{u))R{Lo) , (18) 

where the dot means the time derivative, O is the area (volume in three dimensions) of the domain and 

R{<j>) = R{4>)er , (19) 

with the radial unit vector 6^. The distance from the center of gravity to the interface is denoted by R{4>) with 
the angle (}> with respect to the x axis. We assume that R{(t>) is a single-valued function of (j) for suSiciently weak 
deformations. The infinitesimal length dui along the interface is related with (i0 as 



dR 



d4> 



dcj) = + R'^d(t> , (20) 



where the prime indicates the derivative with respect to (j). The tangential unit vector t and the unit normal n at a 
position on the interface are given, respectively, by 

t = -^^^={R'er + Res) , (21) 

n = ^ {Rp,.-fl'p.A , (22) 

where is the azimuthal unit vector. Prom these expressions, one obtains 

K{ct>) = _ ^ +2^ -J _ ^ , (23) 



VR' + R'' 

RdR/d 
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Deformations of a domain around a circular shape with radius Rq are written as 

R{(l>) = Ro + 5R{(j),t) , (25) 

where 

oo 

5R{4>,t)= Cnity^"^ . (26) 

n= — oo 

Note that since the translational motion of the domain has been incorporated in the variable p, the modes c±i should 
be removed from the expansion (26). The mean curvature and the normal component of the velocity are given up to 
the first order of the deviations, respectively, by 

1 1 °° 

m,t) = -^-^ E (n'-l)cn(^)e*"^ (27) 

^ n— — oo 

OO 

y((/.,i) = vn+ Yl c„(^)e™^ (28) 

n= — oo 

For an isolated domain, the Fourier transform of 6{u — p) = 9{R —\i — p\) is given by 

0q=- I dr exp(-ig • r) . (29) 

Jr<R 

Here the Fourier transformation is defined by 

A{r,t) = [ A,(^)e*'•^ (30) 



Ag{t) = J drA{r,t)e-^'^-\ (31) 

where = J dq/{2iTY with d the dimensionality of space. Substituting (25) into (29) yields up to first order of the 
deviations 

es) = of^ + o^l\ (32) 

where 

= ^.A(5i?o) , 
q 

n 

and J„ [x) is the Bessel function of the n-th order defined by 

/■27r 

/ dee'("^±^™"^' =27ri±"J„(0) . (34) 
Jo 

Here we summarize the formulas for the Bessel function 

J_„(z) = (-l)"J„(z), (35) 

z"J„_i(^) = ±{z^J^{z)) , (36) 
2n 

— Jn{z) = Jr,-l{z) + Jn+l{z) , (37) 
Z 

2J'^{Z) = Jn-l{z) - Jn+l{z) . (38) 

Since the amplitude c„ is complex, it is convenient to introduce the linear combinations of its real and imaginary 
parts. The modes c±2 represents an elliptical shape of domain. We introduce a second rank tensor as follows; 

'S'll = -S'22 = C2 + C_2 , 

'S'12 = S21 = i{c2 - c-2) . (39) 
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For an elliptical domain, R{(j)) is represented as 

i?(</.) =i?o + yCOs2(0-(/.2) , (40) 

where 62 is a positive constant and 02 is the angle between the long axis of the elliptical domain and the x-axis. The 
tensor S can be written in terms of ^2 as Sn = ((^2/2) cos 2(^2 and S12 = ((52/2) sin2(/>2. This is represented in terms 
of the unit normal N = (cos (f)2, sin ^2) along the long axis as 

Sc,p ^ 52{N^Np - ^N^) , (41) 

which is the same as the nematic order parameter tensor in liquid crystals [26] . 

The modes c±3 are necessary to represent the head-tail asymmetry of a propagating domain. If the deformation 
i?((^) is written as 

R{(t)) =Ro + 63 cos 3{(t) -h) , (42) 

with a positive constant ^3, one obtains 

Ti = C3 + C-3 = 63 cos 3(^3 , 

T2 = i{c3 - c-3) = ^3 sin 3(t>3 ■ (43) 



It is convenient to introduce a third-rank tensor associated with the n = ±3 modes as follows; 
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t/a/3.= ^ E iV^ivWiVW, (44) 



m=l,2,3 



where 



ATW = (cos (/)3, sin 03) , (45) 
Ar(2) = ( cos(03 + y), sin(03 + y)) , (46) 
JV(3) = ( cos(03 - y), sin(03 - y)) • (47) 

Prom these definitions, one obtains Um = Ti and f7222 = —T2. Furthermore, it is readily shown that there are 
relations among the components as 

Uni = —U122 = —U212 = —U221 , 

C/222 = -Uu2 = -U121 = -U211 . (48) 

The tensor (44) is the same as the order parameter for banana (tetragonal nematic) liquid crystals in two dimensions 
[41, 42]. 

In the following sections, we shall derive a coupled set of the time-evolution equations for v, c±2 and c±3. However, 
it is impossible to derive the set of equations exactly in a general condition. As mentioned in Introduction, we employ 
an assumption, that is, the smallness of the propagating velocity. Since the circular shape of a domain is assumed 
to be stable when it is motionless, the deformation of the domain is expected to be weak near the drift bifurcation 
where the propagating velocity is small [18]. In this condition, the trancation of the modes up to n = ±3 is justified. 



V. EQUATION OF MOTION FOR p{t) 

The aim of this section is to derive the equation of motion for the center of gravity from the interface equation (13). 
When the velocity is small, one may expand (13) with (14) in powers of V. Up to the third order, one obtains 

2w+ ^ ^ = l-2p. (49) 

The value w of the inhibitor v at the interface is given from Eq. (16) by 

W = W(Q) + U)(i) + W(2) + W(3) , (50) 
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where 

W(0) 



= / G^e^e"^-''^'^^ , (51) 

Jq 

W(l) = W(ll) + W(i2) 

= i [{v q)G2(?,e*'J-«(-) 

Jq 



- |^G^(^)e-«M, (52) 

W(2) = W(21) + W(22) + W^(23) 



Jq 

Jq 



W(3) 

with 



+ y^G^(^)e*'?-«(-), (53) 
= -i [{v q)3G^^ge^<'-«('^) , (54) 

Jq 



Wc have dropped out terms having a factor t;^ in (53) and vvO and the second derivatives in (54). 

In order to obtain the equation for p, we operate (1/il) / dujR{uj) to Eq. (49). As mentioned at the end of the 
preceding section, the basic approximation is the expansion in terms of v and c„ ignoring the higher order time 
derivatives. The validity of these approximations will be discussed in detail in section 7. 

Since the Lagrange multiplier L is independent of the angle 0, one has / dujR{uj)L = J duiR{uj)L^ = 
J duj R{llj)V'^ L = for a circular domain. Furthermore, it will be shown later in section 6 (just after Eq. (91)) 
that L ~ 0{v^), one may ignore / dujR{uj)V L"^ up to the third order of t;. It is also noted that / du:R{uj)K(uj) = 
up to the first order of deformations since the modes c±i are excluded in Eq. (26). The term eVK and the higher 
order terms of e are ignored. As a result, eK + L in the third term on the left hand side of Eq. (49) can be dropped 
out in the derivation of equation for p. 

The zero-th order terms in Eq. (49) which represent a motionless circular domain are given by 

- TT^ + 1 - 2P - 2 / G,0(o) e^'^-^'"' = , (56) 

^^0 Jq 

with = i?(")e,.. The Lagrange multiplier L can be omitted at this order since it is absorbed into the constant p. 
Equation (56) gives us the equilibrium radius Rq of the circular domain as a function, e.g., of the parameter p [39]. 
Since the radius has been fixed by the constraint (6), this implies that the parameter p is not independent but should 
be related to W in (6). 

From the higher order terms in Eq. (49), one obtains in two dimensions 

T 3r^ 

2w + -V - —v\v\^ = Q , (57) 



where 
with 

if(ii) 



2 64 

W = W(ii) + t«(21) + W(3) , (58) 

i|da;|(t;.q)G^^,^e'<'«(-), (59) 



«'(2i) = -T, duj ii,. q)G30,— , (60) 



w 



(3) 



-l/../^(...)G3.,|e-«(" 
-LJdivJiv.qfG%^e^''-^(-K (61) 
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These have been evaluated for a circular domain in two dimensions previously [18]. Substituting those results into 
Eq. (57) gives us 

mi) + i(r - Tc)v + gv\v\'^ = -25w(ix) , (62) 

where 

POO 

m = 2Ro dqqGlJx{qRof , (63) 
Jo 

Tc = 4i?o / dqqGlMqRo f , (64) 
Jo 

9 = ^fdqq'GtMqR^^r-^-^. (65) 

The loft hand side of Eq. (62) has been obtained in rcf. [18]. The right hand side is defined as Sw(^ii) = 'U'(ii) — 

and vanishes identically for a circular domain where if^lil) is given by (60) with 0g = O^q^ and R{uj) = R^°\io). The 
coefficient m is positive and g is shown to be positive for r ~ Tc. Therefore, Eq. (62) indicates a bifurcation (drift 
bifurcation) such that a motionless circular domain = is stable for r > Tc whereas it looses stability for r < Tc 
and undergoes propagation at the velocity = (tc — r)/(2(/). 

The right hand side of Eq. (62) is given up to the order of 0{vCn) by 

+ Lji^j(v. ,)GlC-^e--'"""'> ■ (66) 

The following formula is useful to evaluate Sw(^iiy, 

1 J (i^e'''-^^-) = ^MqRo) 
+ ^ E c„(t)e^"^27ri" J„_i(giio) 

+ ^^c„(t)e*"«27ri"(l-n)J„(gi?o), (67) 

^ n 

where is the angle between q and the x axis, i.e., q = q{cos6,sm6). This formula is valid up to the first order of 

the deformation. 

Now we calculate each term in Eq. (66). The first term is written as 

Sw[ll^^ = Vo^H^^^ , (68) 
where the repeated indices imply the summation and 

h[1^ = -47riio Ec„(Or" / ^Gle'^'^J,{qRo)JniqRo) 

= ai(c2 +c_2) , (69) 



with 



Similarly one obtains 



ai = ^ [ dqq^GlMqRo)J2{qRo) ■ (70) 



Hi^ = -47ri?o ^c„(i)r" 

n 

f M^Gy-'j,{qR^)UqRo) 

J Q Q 



X 

ai{c2-c-2)i. (71) 
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together with the relations = —-^11'' and = . It should be noted that because of the factors in (69) 

and qxqy in (71), only the components n = ±2 contribute to iJ^P and . 
The second term in Eq. (66) can be written by using the formula (67) as 

(2) _ ^ jj(2) 
(ll,/3) ~ ''"-"a/3 



5wfl,, = v^H^^l , (72) 



where 



2 

q"q 



d 



q dqp 



e 



1 — n 

X [qJn-i{qRo) + —^Jn{qB^)] 

= -^y,Cn{t)i" [ e'^'[qJu-i{qRo) 

+ ^^"('^^o)]^|(G^^r)- (73) 

In the derivation of the second line, we have used the fact cq = which comes from the domain area conservation. 
Each component of Eq. (73) is readily obtained as 

H'if = a2(c2 + C_2) , 

= a2(c2-c_2)z , (74) 



where H^f = -H[f> and H[f = H^f and 



1 f°° 

0.2 = / dqq'^[qJi{qR^) 

AttRo Jo 



- ^J2(5i?o)]^(GX°')- (75) 
Putting the expressions of H^J with n = 1 and 2 together, Eq. (62) is finally written as 

+ Tc)Va+gVa\vf = -aV/^S/^a , (76) 

where 

a = 2(ai+a2), (77) 

and the tensor SjSa has been defined in Eq. (41) with (39). The coefficient a is evaluated numerically and is displayed 
as a function of Ro = Rqk in Fig. 1. It is noted that the coefficient a is negative for < Rq < oo. When Rq <C 1, one 
has an analytical expression d = cD'^/Rq = — (3/4)| lni?o|. The numerical result in Fig. 1 is consistent with this. 

VI. EQUATIONS OF MOTION FOR c±2(i) AND c±3(t) 

In order to derive the time-evolution equations for c±2{t) and c±3{t), we start with the interface equation of motion 
(13). Since it will be shown in section 7 that c±2{t) ^ 0{v'^) and c±3(i) ^ 0{v^), the coupling such as c„c„i is 
negligible for slow dynamics of a domain. Therefore, we may consider only the linear order of the deviation 6R. (This 
takes account of the coupling between c„ and p.) Linearizing Eq. (13) with respect to the deformation, one obtains 

d6R e ,d^SR 

where 5w = w — wq with wq the value of w for a circular domain and we have used the fact that Tdc{wo)/dwo = —4 
as shown from Eq. (14). The last term L is necessary in order to eliminate any zero modes which might arise from 
5w. 
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FIG. 1: The coefficient a = aD^ /R(, as a function of Rq = kRq. 



First, we consider the n = ±2 modes. From Eqs. (51), (52) and (53), one notes that there are three terms in 5w, 
which produce n = ±2 modes; 



5w — Swq + Swi + Sw2 , 



where 



J q J q 

6wi — Swii + dwi2 

t j{q ■ v)Gp,e^'' ''^-^ - ^ G^^e"'-^^-) , 

Jq 



Swq ~ 



5W2 = 



(79) 

(80) 

(81) 
(82) 



Here we have ignored the term with d^Oq/dt^ which causes c„ and are higher order in the final equation of motion. 
By using Eq. (33), these are readily evaluated as 



6wo 
6wii 

Swi2 
SW2 



n 

"^(^'l-*^^2)E^»^"e*^"^'^* 

n 

l{vi+iv2)Y.B„c,,e'("-^^^ , 

n 

n 

-^v' + ^[{li-l-v^)-^v^V2) 



(83) 

(84) 
(85) 

(86) 
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where 



/■OO 

A„ = Rq dqq^G'l[Jn{qRQ)Jn+i{qIh) 
Jo 



+ Ji{qRo)g^Ji{qRo)] , (87) 



/■OO 

B„ = i?o / dqq^Gl[-Jn{qRo)Jn-i{qRo) 
Jo 



D. 



+ Ji{qRo)-^Ji{qR^)\ , (88) 

/■OO 

= Ro dqqGg[Ji{qRof - JniqRof] , (89) 
Jo 

/•OO 

En = Ro dqqOlJniqRof , (90) 

JO 

/■OO 

Ge = Ro dqq^GlJt{qRo)Ji+i{qRo) • (91) 
Jo 

Note the relation that A_„ = B„. It is also noted from Eq. (86) that 5w2 contains the n = mode proportional to 
u^, which should be absorbed into the Lagrange multiplier L in Eq. (78). From the n = ±2 modes, one obtains 

Swr^ + 5wr-^ = ^{vl- vl) , (92) 
i(<5<=^ - Sw^=-'^) = GiViV2 . (93) 
After extracting the n = ±2 terms from Eq. (78), equation for the tensor S defined by Eq. (41) is given by 

= -K^S^^ + b[v„V0 - ^-^v^] + biUc,0^v^ , (94) 

where = t — b = — 4Gi, K2 = {3e/Rf^) — 4I?2 and 61 = 2B3. It is noted that the coefficient b is negative. In 
the vicinity of the bifurcation r ~ Tc with Tc given by (64), the coefiicient r2 is positive. The coefficient K2 becomes 
negative for sufficiently large values of Rq indicating an instability of a motionless circular domain [39] . The coefficient 
61 is plotted as a function of Rq in Fig. 2. It is positive for Rq < 0.8 and negative for Rq > 0.8. 

By putting all the n = ±3 terms together, the time-evolution equations for c±3 or Tj (i = 1 and 2) defined by (43) 
take the following form 

= -KsT, 

dt 

+ divi{vf - 3v^) + d2{viSn - V2S12) , (95) 

- diV2{vl - 3t;i) - d2{v2S22 - V1S21) , (96) 

where T3 = t — AE^ and = Se/R"^ — AD^ and the relation = —S22 has been used. The second term in Eqs. 
(95) and (96) arises from W(3-) given by (54) whereas the third term in Eqs. (95) and (96) comes from given by 
(52) and the coeSicients di and d2 are given, respectively, by 

POO 

di = Ro dqq^GpiiqRo)MqRo) , (97) 
Jo 

POO 

d2 = 2Ro / dqq^Gl[j2(,qRo)MqRo) 
Jo 

+ JiiqRo)j^Ji{qRo)] . (98) 
The coefiicient ^2 is displayed as a function of -Rq in Fig- 3. Both di and ^2 are positive. 
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Equations (95) and (96) can be represented in terms of the third-rank tensor Ua^j defined by (44) as 

i 3 — -^3(^a/37 



- -^{^al3S'~ir^ + 5i)^Sari + 5'^aSi3rf)\ . (99) 

Other terms such as U vv are higher order as shown in the next section. 

VII. ORDER ESTIMATION AND STATIONARY SOLUTIONS 

In this section, we shall discuss how to justify the approximations employed in the derivations of the time-evolution 
equations Eqs. (76), (94) and (99). As emphasized, the basic assumption is that we are concerned with the domain 
dynamics near the supercritical drift bifurcation t = Tc- Therefore, the smallness parameter is 6 = Tc — t. It is 
readily shown from Eq. (76) that all the terms are of the order of 5"^/^ since v ~ 0{S^^^) and time is scaled a,s i — tS 
and S ~ 0{6) as can be seen from the first and second terms in Eq. (94). Similarly, one notes from Eq. (99) that 
U^OiS^/^). 

The above order estimation tells us that the terms ignored in Eq. (76) are indeed higher order in 6. For example, 
d^v/dt'^ 0((5^/^) and SSv 0(6^/^). On the other hand, the first and second terms in Eq. (94) are of the order of 
0{S) whereas dS/dt ~ 0{S'^) and Uv ^ 0{S'^) and all other terms ignored arc higher order. Equation (99) has also a 
similar property. Therefore, up to the leading order, the time derivative in S and U should be ignored. This is not 
surprising since, for ^ — > 0, the center of gravity is a slow variable but the deformations around a circular shape are 
not necessarily slow and might be eliminated adiabatically. This fact does not cause any diflaculties in the study of 
the stationary shape of a steadily propagation domain. 

It is more convenient to employ the following representation 

vi = V cos (j), V2 = V sin (f) , (100) 

Sn = ^500826*, 5i2 = ^ssin26' , (101) 

f^lll = C3 + c_3, U222 = -«(c3 - C-3) , (102) 

C3 = |e-3^^3 = f e-3iv^ ^ (103) 

where v should not be confused with the inhibitor v in Eq. (3). Substituting these into Eqs. (76), (94) and (99) yields 

= v{j -v^) -^a'svcos2{e - (j)) , (104) 
j^<j> = -^a'ssm2{0-^) , (105) 

= -KS + 6V cos 2(6» - (j)) , (106) 



dt 2s 



sm2{e-(l)) , (107) 



= — Kz + d[v^ cos 3 {if — (j)) 

+ d'^sv cos(3(/3 -20-^) , (108) 

d d\ r, . „, 

_^ = -^t; sm3(cp-0) 
d' 

- -^susin(3(^-26'-(/>) , (109) 
oz 
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where 7 = (tc — T)/(2m), a' = a/m, k = K2/T2, b' = b/r2, K = d[ = 2di/r3, and d'2 = d2/T^ and we have 

put g/m — 1 without loss of generahty. We have ignored the bi term in Eq. (94) since it is higher order. In what 
follows, we drop the prime in a', 6', d'^ and d'2. 

It is noted that Eqs. (104)-(107) are closed. In the previous paper, we have shown that there are two stable 
stationary solutions depending on the parameters [25]. For simplicity, we restrict ourselves to ab > 0. Here we 
consider a straight motion propagating, without loss of generality, along the x axis, i.e., = 0. When b is positive, 
= 0, and the steady value of the amplitudes is given by w§ = 7/(1 + i3) and sq — bv^j k. where B = ab/{2K). whereas, 
when b is negative, the solution is given hy 6 = 7r/2, = 7/(1 + B) and sq = —bv^/K. 

As shown in the preceding section, the both constants a and b are negative in the reaction diffusion equations studied 
in this paper. Therefore, the long axis of an elliptical domain is perpendicular to the velocity vector. Substituting 
the stationary solutions given above into Eqs. (108) and (109), one obtains the time-independent solution 



Kzq — di^Q cos 3(^0 
= diWQ sinSipo - 



- d2SQVo cos(3i^o - tt) 
d2SQVQ sin(3i^o - tt) . 



Since zq should be positive, we find the following condition irrespective of the sign of 6; 

di + —d2 > =^ cos 3930 = ^ (ySo = , 

K 

di H — c?2 < ^ cos 3(^0 = ^1 => fo ^ , 

K 

and 



di + -d2 

K 



K 



(110) 
(111) 



(112) 
(113) 

(114) 



The stationary shape of a domain propagating along the x axis ((/) = 0) is shown in Fig. 4(a) for tpo = 0. It should 
be noted that in the present reaction diffusion system, the constants di and c?2 are positive whereas b is negative. 
When K is large, the straight motion is stable and the inequality (112) always holds. Therefore, we have the deformed 
domain shown in Fig. 4(a). We emphasize that this is consistent with the numerical simulations in ref. [10]. On the 
other hand, when b is positive (the constant a is also assumed to be positive), the long axis of an elliptical domain 
is parallel to the velocity vector [25]. In this case, the stationary domain shape is displayed in Fig. 4(b) for (^0 = 0. 
If the condition (113) holds, i.e., (p^ ~ tt, the shape of domain is given by the mirror symmetry with respect to the 
vertical axis of those in Fig. 4. 





-Ro Ro 



FIG. 4: (a) The stationary shape of a domain traveling along the x axis (0 — 0) for (po = and (a) 6 = n/2, and (b) 6 — 0. 



Hereafter, we shall discuss the properties of the time-evolution equations (76), (94) and (99) from more general 
point of view. In fact, all the terms are derived by considering the possible couplings which can be generated from 
the first rank tensor (vector) v, the second rank tensor S and the third rank tensor U not necessarily restricted to the 
reaction-diffusion equations. Therefore, if the coefficient r2 ^ 0{S~^) (and T3 ~ 0(1)), one has to retain the dS/dt 
term in r2 whereas the Uv term of the order of 0((5^) can be ignored. The resulting coupled equations (76) and (94) 
are the same as those introduced in ref. [25] where the rotation bifurcation has been found. The bifurcation occurs at 
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such that the stationary straight motion becomes unstable for 7 > 7c. Assuming a circular motion appears for 7 > 7c, 
we put V = Vr, s = Sr, = Lot + and tj) = ojt. Substituting these into Eqs. (104)-(107), one obtains after some 
algebra = 7 — k/2 , = bv^/a, cos( = n/aSr and = {ab/4:){vf — v^) where Vc = K/{abY^'^. It has been shown 
that the frequency lo continuously increases from zero at 7 = 7c. We have carried out the stability of the straight 
motion analytically and that of the circular motion numerically and have found that when the straight motion is 
stable, the circular motion is unstable and vise versa. That is, there is no parameter regime where these two motions 
coexist [25]. 

The stationary shape of a circular motion can also be obtained without any difficulty. There are two types of circular 
motion; clockwise and anti-clockwise rotation. We consider only the anti-clockwise rotation putting = ujt + S/3. 
Then, the set of equations (108) and (109) become 

KZr — diV'^ COS S + d2SrVr C0s{5 — ^) , (116) 
3ZrC0 = —diV^^sinS — d2SrVrSm{S — . (117) 

From these two equations, the stationary solution is given by 

cos^ 

sin 6 

where ds = d2/a and 

tan (5 

From Eqs. (118) and (119), one has 

The replacements u — > —w, ( —(^, and 5 

\ 111. SUMMARY AND DISCUSSION 

In this paper, we have studied the domain dynamics in an excitable reaction-diffusion system with the global 
coupling. The set of time-evolution equations has been derived in the singular limit that the interface width is 
infinitesimal. The main results are Eqs. (76), (94) and (99). The coefficient a in Eq. (76) and b in Eq. (94) are 
negative whereas the coefficients di and d2 in Eq. (99) are positive. Therefore, the deformed shape in the steady 
state takes the form shown in Fig. 4(a)). We emphasize that this result is consistent with the computer simulations 
of Eqs. (2) and (3) with (4) and (5) [10]. 

The present theory is restricted to two dimensions. Extension to three dimensions is possible if one expands 
deformations around a spherical shape in terms of spherical harmonics. Representation of the set of equations by 
introducing some tensor variables as in the present theory is also convenient. However, it is expected that a propagating 
domain in three dimensions are elongated not only to a rod shape or a cone shape (corresponding to the form in Fig. 
4(b)) parallel to the propagating direction but also to a jellyfish-like shape (corresponding to the form in Fig. 4(a)) 
when the elongation is perpendicular to the propagating direction. Therefore, the direct use of the same tensors Sai3 
for the nematic order parameter and Ua/Sj for banana-shape liquid crystals seems not to be possible. 

Apart from the generalization to three dimensions, there are several important open problems. First of all, the 
present theory should be generalized to a collision process of a pair of dcformable self-propelled domains. If the 
interaction is repulsive, the motion of the two domains becomes slow as approaching to each other and hence those 
shapes are expected to change as shown in the present theory. However, the existence of another domain itself tends 
to deform generally the domain, which should be incorporated into our theory. Next, the dynamics of domains 
undergoing the circular motion will be investigated in detail when the global orientational coupling is imposed. 
Since a clock-wise motion and a counter clock-wise motion coexists, the orientational ordering is interrupted so that 
complex dynamics, such as synchronization, localization and chaotic motions appear depending on the parameters 
[25]. Numerical simulations and a theoretical study reducing to non- linear coupled oscillators are now being carried 



{divf 

dfv^ + dzvf. (2Krfi + bd2) ' 
2ijdT,K + 3a; {div"^. + nd^) 
djv^ + davf {2Kdi + 6^2) ' 

2uidiK + 3a; {div'^ + nd^) 
{div"^ + Kda) K - 6uj^d3 ' 

dlv^ + djvf {2Kdi + bd2) 
-5 give us the solution for a clockwise rotation. 



(118) 
(119) 

(120) 
(121) 
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out [43]. Third, the helical motion in three dimensions which is an extension of the circular motion in two dimensions 
must be investigated both numerically and analytically. Finally, the stochastic dynamics of deformable self-propelled 
domains will be developed by adding random noise terms in the equations derived here. As a related study, quite 
recently, Sano at al. [44] have introduced and investigated a dynamical model for the motion of amoeboid cell to 
analyze the experimental results obtained. We hope to publish these investigations somewhere in the near future. 
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